rm(list = ls())
set.seed(1)
library(WhatIf)

Ns=c(30,50,100,200,300,400,500,600,800,1000)
Ks=c(1,2,3,4,5,6,7,8,9,10)
sims=1 
P=.5

N={}
K={}
mean.dist={}
nearest.dist={}
pct.hull={}

for(n in 1:length(Ns)){
  for(k in 1:length(Ks)){
    print("N and K:")
    print(c(Ns[n],Ks[k]))
    pct.hull.t <- {}
    mean.dist.t <- {}
    nearest.dist.t <- {}
    for(s in 1:sims){
      print("sim:")
      print(s)
      x <- matrix(rnorm(Ns[n]*Ks[k]), Ns[n], Ks[k])
      treat <- as.numeric(runif(Ns[n])>P)
      x0 <- x[treat==0,]
      x1 <- x[treat==1,]
      if(Ks[k]==1){
        x0 <- as.matrix(x0)
        x1 <- as.matrix(x1)
      }
      h0 <- whatif(data=x0, cfact=x1, distance="euclidian",
                   return.distance=T)
      h1 <- whatif(data=x1, cfact=x0, distance="euclidian",
                   return.distance=T)
      pct.hull.t[s] <- (sum(h0$in.hull)+sum(h1$in.hull))/(nrow(x0)+nrow(x1))
      nearest.dist.t[s]=mean(rbind(as.matrix(apply(sqrt(h1$dist),1,min)),
                      as.matrix(apply(sqrt(h0$dist),1,min))))
      h <- whatif(data=x, cfact=x, choice="distance", distance="euclidian",
            return.distance=T)
      mean.dist.t[s]=mean(sqrt(h$dist))
    }
    N=rbind(N,Ns[n])
    K=rbind(K,Ks[k])
    mean.dist=rbind(mean.dist,mean(mean.dist.t))
    nearest.dist=rbind(nearest.dist,mean(nearest.dist.t))
    pct.hull=rbind(pct.hull, mean(pct.hull.t))
  }
}

print("N,K,mean.dist,nearest.dist,pct.hull")
Dat=cbind(N,K,mean.dist,nearest.dist,pct.hull)

source("compact.R")
sink("NK.out", append=T)
print(compact(Dat))
sink()
detach(package:WhatIf)
detach(package:lpSolve)
gc()
q()
